#####################################################################
# Code used to create example data included with `ResistanceGA`     #
# Written by Bill Peterman                                          #
# 24 June 2014                                                      #
#####################################################################

library(RandomFields)
library(ResistanceGA)
rm(list = ls())
set.seed(123)

# Random point and raster parameters
r.dim <- 50       # number of cells on a side
cell.size <- .5        # raster cell dimension
min.point <- floor(0.25*(r.dim*cell.size))#+ (cell.size/2)       # minimum coordinate for generating random points (mult by 0.25 to prevent edge effects)
max.point <- (r.dim*cell.size)-min.point        # maximum coordinate for generating random points

# Number of "Sample locations" to generate. This example will generate points on a square grid, so choose a number that has an even square root
n <- 25
x <- ceiling(seq(from=min.point,max.point, length.out=5)) + (cell.size/2)     # set x & y range to draw random samples from
y <- floor(seq(from=min.point,max.point, length.out=5)) + (cell.size/2)
Sample.points<-expand.grid(x,y)

Sample.coord <- SpatialPoints(Sample.points)
samples <- cbind((1:n),Sample.coord@coords)
colnames(samples)<-c("pop","x","y")

# Specify random model
model<-RMexp() +
  RMtrend(mean=10)

grid.vars <- GridTopology(cellcentre.offset=c(cell.size/2, cell.size/2),
                          cellsize=c(cell.size, cell.size),
                          cells.dim=rep(r.dim,2))

rf.sim <- RFsimulate(model, x=grid.vars,n=2) # Simulate 2 surfaces

# The first layer is continuous
continuous <- raster(rf.sim[1]) # Define the first as a continuous surface
names(continuous)<-"continuous"
plot(continuous)
plot(Sample.coord, pch=16, col="blue", add=TRUE)

# The second layer will be converted into a 3-class categorical surface
categorical <- raster(rf.sim[2])
names(categorical) <- "categorical"
cat.cut <- summary(categorical) # Define quartiles, use these to define categories

categorical[categorical<=cat.cut[2]] <- 0
categorical[categorical>0 & categorical<=cat.cut[4]] <- 1
categorical[!categorical%in%c(0,1)] <- 2
plot(categorical)
plot(Sample.coord, pch=16, col="blue", add=TRUE)

# Create linear feature class
feature <- matrix(0,r.dim,r.dim)
feature[27,] <- 1
feature[,22:24] <- 1
feature <- raster(feature)
extent(feature)<-extent(categorical)
names(feature)<-"feature"
plot(feature)
plot(Sample.coord, pch=16, col="blue", add=TRUE)

# Make a RasterStack
resistance_surfaces <- stack(categorical,continuous,feature)